16th Computational Fluid Dynamics Conference, 23-26 June 2003, Orlando, FL 


Failure of Anisotropic Unstructured Mesh 
Adaption Based on Multidimensional Residual 

Minimization 

William A. Wood* and William L. Ivleb* 

NASA Langley Research Center, Hampton, VA 23681 

An automated anisotropic unstructured mesh adaptation strategy is proposed, imple- 
mented, and assessed for the discretization of viscous flows. The adaption criteria is based 
upon the minimization of the residual fluctuations of a multidimensional upwind viscous 
flow solver. For scalar advection, this adaption strategy has been shown to use fewer grid 
points than gradient based adaption, naturally aligning mesh edges with discontinuities 
and characteristic lines [Wood, AIAA 99-3254]. The adaption utilizes a compact stencil 
and is local in scope, with four fundamental operations: point insertion, point deletion, 
edge swapping, and nodal displacement. Evaluation of the solution-adaptive strategy is 
performed for a two-dimensional blunt body laminar wind tunnel case at Mach 10. The 
results demonstrate that the strategy suffers from a lack of robustness, particularly with 
regard to alignment of the bow shock in the vicinity of the stagnation streamline. In 
general, constraining the adaption to such a degree as to maintain robustness results in 
negligible improvement to the solution. Because the present method fails to consistently 
or significantly improve the flow solution, it is rejected in favor of simple uniform mesh 
refinement. 


Nomenclature 

A, A Flux Jacobian in conserved or auxiliary variables 

a Speed of sound 

F Flux of conserved quantities 

/ Identity matrix 

i, j Cartesian unit vectors 

£ Edge length 

M Upwinding matrix 

h Outward-normal unit vector 

P Pressure 

S Area 

s Entropy: ds = dp — dP /a 2 

T Temperature 

U Conserved variables 

u,v Velocity components 

V Velocity vector 

W Auxiliary variables: dW = (ds,pdV,dP) 
x, y Cartesian coordinates 

Z Parameter vector: ^/p(l,u,v,H) 

a, p Generalized wave speeds in directions 

7 Ratio of specific heats 

£ Curvilinear coordinate on triangle edge 12 

y Curvilinear coordinate on triangle edge 23 

T Minimization functional 

<f>, <j) Fluctuation in conserved or auxiliary variables 
p Density 

5 Functional weighting matrix 

fl Integration element 
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Superscripts: 

( ) T Transpose 

( ) ?! Viscous component 

( ) x,y Spatial component of a vector 

Subscripts: 

( ),: Node number 

( )t Triangle 

Over-bars indicate linear averages and tildes indi- 
cate average values assuming linear variation of the 
parameter vector. 

Variables lacking explicit dimensions are in nondi- 
mensional form, normalized by the free stream values 
of density, velocity, viscosity, and temperature along 
with a characteristic length, as appropriate. 

Introduction 

T HOMAS et al . 1 identify viscous flow unstructured 
mesh adaptation as a critical enabling technology 
for the inclusion of high fidelity computational fluid 
dynamics in the vehicle design cycle. The development 
of robust automated unstructured grid generation and 
adaption methods will directly reduce the time re- 
quired for both the initial grid generation for complex 
vehicles and the solution computation for complex flow 
fields. Notable achievements in this field of unstruc- 
tured mesh adaption have recently been demonstrated 
by Venditti and Darmofal , 2 using adjoints to optimize 
the mesh distribution, and by Habashi et al ., 3,4 using- 
anisotropic feature adaption. 

Concurrent work in flow solvers for hypersonic 
aerothermodynamic applications has shown promise 
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for nonlinear multidimensional upwind fluctuation 
splitting schemes as being more accurate and less 
grid sensitive than traditional locally one-dimensional 
flux splitting schemes. 5-7 Unstructured mesh adap- 
tion based upon the minimization of the fluctuation 
splitting distribution has shown some remarkable nat- 
ural alignment properties for linear 8,9 and nonlinear 10 
model problems, but has not been demonstrated for 
the full system of equations for viscous fluid dynam- 
ics. Alignment of the mesh to the bow shock has been 
shown to be a critical factor in the accuracy of hy- 
personic solutions, 11 and the successful extension of 
the natural alignment tendency inherent in fluctuation 
minimization to a high speed fluid dynamics capability 
could be the desired breakthrough to the cost barrier 
of including computational fluid dynamics in the vehi- 
cle design cycle. 

This paper presents an automated anisotropic un- 
structured mesh adaption process based upon the 
minimization of distributions for a nonlinear multidi- 
mensional upwind fluctuation splitting scheme. The 
adaption is local in scope and thus convenient to par- 
allelize. The demonstration case is for perfect gas 
two-dimensional Mach-10 flow over a sting-mounted 
entry capsule profile. 

Flow Solver 

The flow solver is a nonlinear monotone second- 
order accurate multidimensional upwind fluctuation 
splitting scheme as introduced by Sidilkover 12 with ex- 
tensions to viscous flows. 13 The solver produces exact 
convective solutions along characteristics on unstruc- 
tured grids where one edge of each cell is aligned with 
the cell-averaged transport direction. 10 This character- 
istic alignment property of the flow solver is the basis 
for driving the mesh adaption strategy described in the 
following section. Traditional locally one-dimensional 
flux splitting solvers do not share this method-of- 
characteristics property. 

While a full description of the second-order accu- 
rate viscous flow solver is available in Ref. 6, including 
verification and validation results, the first-order dis- 
cretization is briefly sketched here to provide context 
for the residual derivatives that follow in the adapta- 
tion section. The equations of motion are written in 
vector form as 


U t = -V-E + V-F u (1) 

where F and F v are the inviscid and viscous fluxes 
and U is the vector of conserved variables. Integrat- 
ing over a triangular element, Figure 1, the right-hand 
side (RHS) can be expressed in terms of inviscid and 
viscous fluctuations, 0 and <j> v , as 

RHS = - / S7-Fdtt + / V-F v dfl = <t> + (p v (2) 
Jn Jn 



Fig. 1 Elemental triangular domain for fluctuation 
splitting. 


The viscous fluctuation is assembled from contribu- 
tions associated with each node of the triangle, 

JL f . ~ 

4> v =^2<t>i and <Pi=^F v -fii (3) 

i= 1 


The edge of the triangle opposite to node i has length 
l{ and outward unit normal h,;. 

Expressing the inviscid flux in terms of its Jacobian, 
dF = A dU , the inviscid fluctuation can be decom- 
posed along two edges of a triangle, denoted £ for 12 
and for 23, as 

4 > =0 i +v’ 1 (4) 

ft = -fill'll- A (ffl = ^i 3 n 3 -AA r] U (5) 


A transformation to auxiliary variables, dU = U\y dW, 
simplifies the expression for the Jacobian, 
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Upwinding is achieved through the introduction of 
an artificial dissipation fluctuation, 


d>' S = sign(ni-A)0 t = sign( — ng (11) 


which is linearly combined with the basic fluctuation 
to complete the inviscid contribution to the nodal up- 
dates. 
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Adaption Strategy 

The present adaption strategy for two-dimensional 
unstructured meshes with triangular elements is local 
in scope and allows anisotropic stretching. Adaption 
cycles are performed as sequential sweeps over the 
domain performing the basic operations: node dele- 
tion, edge swapping, nodal displacement, and node 
insertion. The process is automated and a variety of 
stopping criteria can be considered, such as a maxi- 
mum grid size, total number of cycles, or minimum to 
maximum residual ratio. 

A node is marked for deletion if the fluctuations, 
both in viscid and viscous , a in all surrounding cells are 
below a threshold. The threshold is customizable, and 
a typical choice would be 2-3 orders of magnitude be- 
low the average nodal value. A trial reconnection of 
the mesh is proposed without the node and the ele- 
ment fluctuations are recomputed, from the existing 
solution, to verify that the fluctuations remain small. 
If so, the node deletion is finalized. 

A node is added at the midpoint of an edge, creating 
two more cells and three more edges, if the fluctuations 
in the cells to either side of the edge exceed a threshold. 
When comparing cell fluctuations, the L 2 -norm of the 
fluctuations from each of the governing equations are 
used, weighted by the inverse of the square root of the 
cell area. 

Edges are flagged to be swapped when the root mean 
square (RMS) of the fluctuations in the cells to either 
side are excessive. If the swapped edge maintains a 
valid grid then the RMS of the fluctuations is recom- 
puted, and the swapped configuration is retained if the 
RMS has decreased. 

Nodal displacements are driven by the minimiza- 
tion of a functional that can be expressed in terms of 
either the conserved or auxiliary fluctuations, as the 
minimization of one implies the minimization of the 
other. Because the goal of the nodal displacements is 
to achieve alignment of the mesh with inviscid disconti- 
nuities and characteristic lines, the viscous fluctuation 
is omitted from the displacement functional. 

The nodal displacement functional is constructed for 
each node i as 


( 12 ) 

T 


where the summation is over all triangles connected at 
node i and S is the weighting matrix. The derivatives 


a The sum of the absolute values of the viscous distributions 
is used, as the net cellular viscous fluctuation is identically zero. 


of X are formed using the chain rule as 
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Having defined the gradient of the functional, the 
method of steepest descent can be applied directly to 
drive the nodal displacements. 

The weighting factor Sy is a symmetric positive- 
definite matrix. Weighting each equation equally with- 
out regard to cell sizes results in Sy = I, while an 
inverse area weighting yields Sy = ^ -I. The deriva- 
tives of Sy for these and other area-weighted choices 
are available in Ref. 14. The present study is con- 
ducted using equal weighting. 

The discrete cell fluctuation, Eqn. 7, can be equiva- 
lently written as 

1 3 

4> (15) 

i= 1 


where Wj = WgZj. The summation is over the three 
nodes defining the triangle. 

The derivatives of the fluctuation with respect to 
the movement of only one node, for example node 2 of 
the triangle, are 


d(p 1 

dx -2 2 
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The derivatives of the fluctuation with respect to 
movements of node 1 or 3 of the triangle follow di- 
rectly as cyclic permutations. 
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The flux Jacobian in auxiliary variables, Eqn. 10, 
has the approximate derivative 


oa ~ av 

~ ~dx' • 


(18) 


and similarly for the y derivatives. It is assumed that 
moving a node does not change the solution at the 
other two nodes of the triangle, so that the variation 
of the cell-average flux Jacobian scales like one-third 
the variation of the velocity at the node being moved, 


dA 

dx 2 


l dV, 
3 8x2 


I. 


(19) 


Similarly, the variation of the Jacobian of the trans- 
formation scales like one-third of the nodal variation, 
lyT ~ I ffj , and is neglected as a sub-principle term 
relative to the change in solution value, 


dWj dZj 

~ Wr ? - 


0X2 


8x2 


( 20 ) 


Since the solution is locally assumed to vary only at 
the node being moved, i.e., 4/0 = yr0 = 0 

x dx 


E^ - w/ y/l 


j =i 


dx 2 


( 21 ) 


The remaining term to evaluate is jZs- . The steady- 
state distribution of the residual fluctuation to the 
node can be written 

E + M q ) a - (I + M /3 ) 0] W 2 Z 2 

T 

= E [< 7 + M «) aW z z i -(! + M /3 ) 0WgZ 3 } , 

( 22 ) 


or 


E K - P + ) WzZ 2 = E {a + w 2 z 1 - p+w 2 z 3 ) , 

T T 

(23) 
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(24) 

(25) 


The upwinding matrices M a and M /3 expressed in 
auxiliary variables reduce to the identity matrix for 
supersonic flow, see Ref. 13. Differentiating Eqn. 23 
while freezing the Jacobians leads to 
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A similar expression can be formed for the y derivative. 

Neglecting the subsonic blending on M a and M /3 
allows a+ and 0+ to be expressed as 


a+ = 


I a, V a >0 fj+ = f 9, % > 0 
1 0, V a <0 ’ p 1 o, %<0 


(29) 

where V a = V -hi and V /3 = V -n 3 . The derivatives 
follow directly, 
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Having determined P- 2 follows from 
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This completes the derivation of all the terms re- 
quired to evaluate the derivatives of the two-dimen- 
sional objective function, Eqns. 13 and 14. Reference 6 
presents the analogous derivation in axisymmetric co- 
ordinates. 
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Fig. 2 Computational domain for Mars Pathfinder 
capsule. 

Results 

Case Definition 

Evaluation of the adaption strategy for hypersonic 
flow is simplified by choosing a two-dimensional lami- 
nar perfect gas case. The test conditions and geometry 
correspond to the Mach-10 wind tunnel tests of Hol- 
lis. 15 The sting-mounted capsule geometry, depicted 
in Figure 2, is modeled as two-dimensional instead 
of axisymmetric for the present study. b The orig- 
inal axisymmetric wind tunnel model of the Mars 
Pathfinder capsule consists of a spherically-blunted 70- 
degree sphere-cone. The body radius is 1 in. and the 
nose radius is jj in. The shoulder radius is in. and 
the aft body angle is 40 degrees. The base radius is 
| in. The sting is 4 in. long with a radius of || in. 

The free stream conditions for the NASA Lang- 
ley 31-Inch Mach 10 Air Tunnel, corresponding to 
a nominal Reynolds number per foot of 0.5 x 10 6 , 
are: P = 69 Pa, T = 53 K, p = 0.0045 kg/m 3 , and 
u = 1416 m/s. The wall temperature is taken to be 
a uniform 300 K. Historical experience with this tun- 
nel at these conditions indicates that laminar perfect 
gas calculations are adequate for comparison with the 
experimental data. 

The metric for accuracy is chosen to be surface heat 
transfer rates. The benchmark solution is obtained 
using the cell-centered flux difference splitting second- 
order accurate Laura code of Gnoffo. 16 A structured 
grid with 125 stream-wise nodes and 513 nodes from 
the surface to the outer boundary, for a total of 64,125 
nodes, is used to generate the benchmark solution. 
The grid is adapted for alignment with the bow shock 
and clustered to the boundary layer. The benchmark 
solution is shown to be grid converged on this mesh, 
aside from small changes at the stagnation point, in 

b Reference 6 reports the axisymmetric results for this case, 
with similar observed behavior and conclusions as for the two- 
dimensional results, but with the added complexity of special 
treatment for the axis singularity. 


Ref. 6. The same reference also shows the grid conver- 
gence of the present solver for this case using uniform 
mesh refinement, converging in agreement with the 
benchmark solution. 

Adaption Component Evaluation 

The individual adaption operations, deletion, swap- 
ping, moving, and insertion, are first tested separately 
to gage their relative merits. Each evaluation begins 
from a converged solution on triangulations of coars- 
ened versions of the structured Laura grid. After 
adaption the solution is re-converged. The adaptation 
will be considered successful if the resultant surface 
heat transfer rates are closer to the benchmark data. 

For point deletion the starting mesh is a triangu- 
lation of 125 x 257 structured grid, containing 32,125 
nodes. The objective for successful point deletion is 
to remove points without altering the solution. For 
the first pass, the adaption is set to remove 10 per- 
cent (3284) of the nodes, producing minimal change 
to the solution. The nodes are predominantly re- 
moved from the free stream. Then a further 10 percent 
(3157) of the nodes are removed, and the heating 
rates remain essentially unchanged, Figure 3, but the 
fluctuation splitting scheme has trouble maintaining 
the bow shock capture due to the further loss of free 
stream points, as shown in Figure 4. The removal of 
more points causes rapid solution deterioration as the 
scheme fails to properly capture the bow shock on the 
coarsened meshes. 

For evaluating edge swapping the starting mesh is 
formed from a triangulated 125 x 129 structured grid 
with 16,125 nodes and 47,868 edges. Edge swapping 
with the fluctuation minimization strategy is not ben- 
eficial for this case. Only about 1 percent of the edges 
are swapped, primarily at the shock, but the result- 
ing solution does not converge due to ringing of the 
bow shock near the stagnation point. Figure 5 shows 
pressure contours in the stagnation point region over- 
laid upon the swapped mesh. Notice the irregular 
contours downstream of the shock. The difficulty in 
this solution is that the bow shock is making discrete 
jumps across grid lines, induced by swapped edges at 
the shock. The unadapted bow shock more closely 
follows the grid lines at the y = 0 symmetry plane. 
Recall that this mesh was originally well-aligned with 
the bow shock from the structured-grid solver, and the 
effect of edge swapping on a more random unstruc- 
tured mesh may be different. Also, on a finer mesh 
the smaller grid spacing at the shock might lessen the 
detrimental effects of discrete jumps in the shock loca- 
tion. The edge swapping assessment is tried again with 
more conservative thresholds, this time only swapping 
half a percent of the edges, but still produces the same 
disappointing results. 

Nodal displacements are evaluated using the 16,125- 
node mesh. The adaption moves 1827 nodes a total 
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a) Heat shield. 



b) Aft body and sting. 

Fig. 3 Surface heating after coarsening, 

solid=unadapted, dashed=10% removed, dot- 
ted=20% removed. 



Fig. 4 Loss of shock capture due to coarsening, 
pressure contours. 



- 0.5 0 


X, in. 

Fig. 5 Shock bulging at stagnation point after edge 
swaps, pressure contours and grid. 

distance of 0.7488 in., for an average of 0.0004 in. and 
RMS of 0.0015 in. This level of movement represents 
a perturbation to the mesh, with the average move- 
ment being an order of magnitude smaller than the 
edge lengths at the shock. The mesh movement occurs 
at the shock, near the shoulder, and in the fore body 
boundary layer. Surface heating rates for this case are 
shown in Figure 6. The solution is worsened on the 
fore body while not much change is seen in the wake. 
Investigating the solution reveals that a wind side vor- 
tex pattern has emerged, Figure 7, that is causing the 
unexpected heat transfer rates on the heat shield. The 
solution is run a further 800,000 iterations with no 
change to this wind side vortex pattern. Note that the 
solver includes eigenvalue limiting, which is known to 
suppress the carbuncle phenomenon for the structured 
grid Roe schemes. The wind side recirculation is not 
a physically accurate prediction nor does it appear to 
match the classic carbuncle pathology, but is rather 
an artifact of discrete mesh point jumps in the bow 
shock location, seen in Figure 7 as choppiness in the 
bow shock contour lines. For this case the adaption 
perturbs a mesh that was initially well aligned to the 
bow shock in such a way as to largely destroy that 
alignment, to the detriment of the solution. 

Point insertion is tested on the unadapted 16,125- 
node mesh. The adaptation adds 3299 nodes, approx- 
imately a 20 percent increase, with heating results 
shown in Figure 8. Many of these new points are added 
to the wind side boundary layer. On the heat shield 
the heating levels rise markedly toward the benchmark 
result, but exhibit a high-frequency oscillation of sig- 
nificant amplitude. There is little change in the wake 
heating. 

Full Adaption Cycle Evaluation 

For the full adaption the starting solution is taken 
on the triangulated 125 x 129 mesh with 16,125 nodes. 
The intention is to look for an improvement in the 
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a) Heat shield. 



b) Aft body and sting. 

Fig. 6 Surface heating after nodal displacements, 
solid=unadapted, dashed=10% moved, dash-dot- 
dot =benchmark. 



Fig. 7 Wind Side vortices produced by nodal dis- 
placements, streamlines and pressure contours. 



s, in. 

a) Heat shield. 



b) Aft body and sting. 

Fig. 8 Surface heating after point insertion, 
solid=unadapted, dashed=20% inserted, dash-dot- 
dot =benchmark. 

coarse mesh solution without increasing the number of 
nodes. Although not all of the component adaptation 
steps were seen to be beneficial, a complete cycle will 
include each step to allow for synergism between the 
components. The strategy for an adaption cycle is to 
delete 10 percent of the nodes, swap 1 percent of the 
edges, move 5 percent of the nodes, and then insert 
back in 10 percent of the nodes. The solution is re- 
converged between each step of the adaption cycle. 

The fluctuation minimization adaption successfully 
removes 1601 nodes, but runs into trouble again while 
swapping. Wind side vortices are spawned in the stag- 
nation region by an oscillating bow shock. In an effort 
to damp the solution the CFL number is reduced by 
an order of magnitude, without producing an improve- 
ment in the solution or eliminating the vortices. The 
number of edges to be swapped is then reduced to 
4 percent, but the same vortices and oscillating bow 
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s, in. 

a) Heat shield. 



b) Aft body and sting. 

Fig. 9 Surface heating after full adaption cy- 
cle, solid=unadapted, dashed=adapted, dash-dot- 
dot =benchmark. 

shock appear. To proceed with this case the edge 
swapping is omitted entirely and the adaption cycle 
continued with the nodal displacement step, where 
706 nodes are moved a total distance of 0.7 in. Fi- 
nally, 1638 nodes are added, yielding the results of 
Figure 9. The fore body heating is generally im- 
proved toward the benchmark solution, although there 
is a high-frequency oscillation in the data starting at 
s = 0.6. The aft body heating is improved to match 
the benchmark between s = 1.2-1. 8, although the orig- 
inal solution was not far off from the benchmark in this 
region. The sting heating is only slightly changed aft 
of s = 3. 

Concluding Remarks 

An anisotropic unstructured mesh adaption strategy 
for hypersonic viscous fluid flows is presented. The 
strategy is implemented for two-dimensional perfect 


gas laminar physics as a prototype for the purpose of 
assessing the benefits of the scheme prior to pursuing 
in a three-dimensional reacting gas implementation. 
The adaption criteria is based upon the minimization 
of residual fluctuations defined by a multidimensional 
upwind distribution scheme for the flow solver. Par- 
ticular emphasis is placed on deriving the derivatives 
of the nodal displacement objective function for the 
upwind fluctuation splitting distribution algorithm. 

The hypersonic demonstration case is a two-dimen- 
sional sting-mounted capsule at laminar Mach 10 wind 
tunnel conditions, modeled with perfect gas air. Sur- 
face heat transfer is the quantity of interest used to 
measure the success of the solution adaptive process. 

Four elementary operations are performed during 
the adaption: point deletion, edge swapping, nodal 
displacement, and point insertion. The limits of ef- 
fectiveness for each of the four operations is probed 
independently to guide their relative weighting in the 
full adaption strategy. Coarsening of the mesh works 
well removing up to about 20 percent of the nodes, but 
leads to a loss of the bow shock capture and a rapidly 
deteriorating solution for more aggressive coarsening. 
The other three adaption techniques generally have 
minimal or negative impacts on the solution, with 
the exception of point insertion where solution im- 
provement occurs but at the cost of high frequency 
oscillations. Overall a lack of robustness is demon- 
strated, usually caused by distortions to the bow shock 
near the stagnation point. Common failure modes are 
an oscillating shock that sheds wind side vortices, or 
a steady shock kink that produces unexpected pres- 
sure contours and surface heating spikes. For the full 
adaption cycle the fore body heating improves but be- 
comes oscillatory, while the aft body and sting heating 
predictions improve slightly. 

One likely cause of the lack of robustness exhibited 
by the adaption scheme at the stagnation point bow 
shock is the conflict between characteristic directions 
between the supersonic upstream flow, which is per- 
fectly upwinded by the flow solver, and the subsonic 
shock layer, where the acoustic characteristics approx- 
imated as discrete waves. A second contributor to the 
lack of robustness is the local scope of the adaption 
strategy. The scheme is intentionally local to be syner- 
gistic with the flow solver which has compact support, 
sharing the same software data structures and loop 
controls. The local scope also is desirable for paral- 
lel processing. However, a larger stencil, an implicit 
scheme, or an adjoint equation criteria may be neces- 
sary to provide sufficient smoothness at shocks. 

Although heat transfer is a variable of primary in- 
terest for aerothermodynamics and may be a first 
consideration for driving the adaption criteria, the 
present results emphasize that for blunt-body hyper- 
sonic aerothermodynamics both the bow shock loca- 
tion and smoothness are critical to the computation of 
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the fore body heat transfer rates. Due to the strength 
of the hypersonic bow shock, discrete jumps in shock 
position from one node to another produce significant 
pressure disturbance waves and entropy variations that 
strongly affect the flow field in the subsonic bubble and 
overwhelm more subtle adaption within the boundary 
layer. 

To summarize, the present adaption techniques do 
not significantly improve the solutions, the adaption 
effectiveness is not consistent, and the solutions dis- 
play a lack of robustness. Point deletion works the 
best, but is risky if applied too aggressively. Uniform 
refinement is both simpler and more reliable than the 
solution adaptive strategy embodied by the present 
method. 
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